#!/bin/bash
set -e

function info() {
echo Usage: `basename $0` in.vcf
exit 65
}

while getopts  ":l:p:" opts
do
        case  $opts  in
        l) bed=$OPTARG;;
		p) out_prefix=$OPTARG;;
		*) info;;
        esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi


. $var

nf=`cat $1|grep '^#C'|awk '{print NF}'`

if test $nf -gt 18; then inb='-an InbreedingCoeff'; fi

echo;echo;echo gatk VariantRecalibrator snp
java $tmp -Xmx$java_memory -jar $gatk \
-T VariantRecalibrator \
-R $ref_genome \
-input $1 \
-recalFile $out_prefix.vcf.recal.vcf \
-tranchesFile $out_prefix.vcf.tranches.txt \
-rscriptFile $out_prefix.vcf.recal.r \
-mode SNP \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 $data_path/gatk/vcf/hapmap_3.3.${genome_assembly}.vcf \
-resource:omni,known=false,training=true,truth=true,prior=12.0 $data_path/gatk/vcf/1000G_omni2.5.${genome_assembly}.vcf \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 $data_path/gatk/vcf/1000G_phase1.snps.high_confidence.${genome_assembly}.vcf \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $data_path/ncbi/dbsnp/${genome_assembly}/All_20150605.vcf \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \
$inb -an QD \
-an FS \
-an MQ \
-an MQRankSum \
-an ReadPosRankSum \
-an SOR \
--maxGaussians 4

# -an DP \
# -U ALLOW_SEQ_DICT_INCOMPATIBILITY #
# -an InbreedingCoeff #population 

# 3 points
#1 InbreedingCoeff #population # need at least 10 samples
#2 exomes
	#2a Depth of coverage no need?
	#2b 30 samples at least # gvcfs?
	#2c --maxGaussians 4

# Depth of coverage
#should not be used when working with exome datasets since there is extreme variation in the depth to which targets are captured! In whole genome experiments this variation is indicative of error but that is not the case in capture experiments.

# InbreedingCoeff
# The InbreedingCoeff is a population level statistic that requires at least 10 samples in order to be computed. For projects with fewer samples, or that includes many closely related samples (such as a family) please omit this annotation from the command line.




echo;echo;echo gatk ApplyRecalibration snp
java $tmp -Xmx$java_memory -jar $gatk \
-T ApplyRecalibration \
-R $ref_genome \
-input $1 \
-recalFile $out_prefix.vcf.recal.vcf \
-tranchesFile $out_prefix.vcf.tranches.txt \
-o $out_prefix.snp_recal.vcf \
-mode SNP \
--ts_filter_level 99.5



echo;echo;echo gatk VariantRecalibrator indel
java $tmp -Xmx$java_memory -jar $gatk \
-T VariantRecalibrator \
-R $ref_genome \
-input $out_prefix.snp_recal.vcf \
-recalFile $out_prefix.snp_recal.vcf.recal.vcf \
-tranchesFile $out_prefix.snp_recal.vcf.tranches.txt \
-rscriptFile $out_prefix.snp_recal.vcf.recal.r \
-resource:mills,known=true,training=true,truth=true,prior=12.0 $data_path/gatk/vcf/Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
$inb -an QD \
-an FS \
-an SOR \
-an MQRankSum \
-an ReadPosRankSum \
-mode INDEL \
-tranche 100.0 -tranche 99.9 -tranche 99.0 -tranche 90.0 \
--maxGaussians 4

# -an DP \
# -an InbreedingCoeff \

echo;echo;echo gatk ApplyRecalibration indel
java $tmp -Xmx$java_memory -jar $gatk \
-T ApplyRecalibration \
-R $ref_genome \
-input $out_prefix.snp_recal.vcf \
-mode INDEL \
--ts_filter_level 99.0 \
-recalFile $out_prefix.snp_recal.vcf.recal.vcf \
-tranchesFile $out_prefix.snp_recal.vcf.tranches.txt \
-o $out_prefix.recal.vcf

 
. $cmd_done






# SNPs and indels must be recalibrated in separate runs (but it is not necessary to separate them into different files). Mixed records are treated as indels.
# The values used in the example above are only meant to show how the command lines are composed. They are not meant to be taken as specific recommendations of values to use in your own work, and they may be different from the values cited elsewhere in our documentation. For the latest and greatest recommendations on how to set parameter values for you own analyses, please read the Best Practices section of the documentation, especially the FAQ document on VQSR parameters.
# Whole genomes and exomes take slightly different parameters, so make sure you adapt your commands accordingly! See the documents linked above for details.
# If you work with small datasets (e.g. targeted capture experiments or small number of exomes), you will run into problems. Read the docs linked above for advice on how to deal with those issues.
# In order to create the model reporting plots Rscript needs to be in your environment PATH (this is the scripting version of R, not the interactive version). See http://www.r-project.org for more info on how to download and install R.

# The InbreedingCoeff is a population level statistic that requires at least 10 samples in order to be computed. For projects with fewer samples, or that includes many closely related samples (such as a family) please omit this annotation from the command line.

# Depth of coverage (the DP annotation invoked by Coverage) should not be used when working with exome datasets since there is extreme variation in the depth to which targets are captured! In whole genome experiments this variation is indicative of error but that is not the case in capture experiments.


# In our testing we've found that in order to achieve the best exome results one needs to use an exome SNP and/or indel callset with at least 30 samples. For users with experiments containing fewer exome samples there are several options to explore:

# Add additional samples for variant calling, either by sequencing additional samples or using publicly available exome bams from the 1000 Genomes Project (this option is used by the Broad exome production pipeline). Be aware that you cannot simply add VCFs from the 1000 Genomes Project. You must either call variants from the original BAMs jointly with your own samples, or (better) use the reference model workflow to generate GVCFs from the original BAMs, and perform joint genotyping on those GVCFs along with your own samples' GVCFs with GenotypeGVCFs.

# You can also try using the VQSR with the smaller variant callset, but experiment with argument settings (try adding --maxGaussians 4 to your command line, for example). You should only do this if you are working with a non-model organism for which there are no available genomes or exomes that you can use to supplement your own cohort.